Automatic detection of rest/activity periods with pyActigraphy

The pyActigraphy package implements several rest/activity detection algorithms.

Broadly speaking, they can be divided in two categories:

  • Epoch-by-epoch rest/activity scoring (inspired by PSG):

    • Cole-Kripke’s, Sadeh’s, Scripps’ and Oakley’s algorithms

  • Detection of consolidated periods of similar activity patterns

    • Crespo’s, Roenneberg’s algorithms

All algorithms have been implemented to return a binary time serie (0 being rest or activity depending on the definition made in the original article…)

In [1]:
import pyActigraphy
/usr/local/lib/python3.7/site-packages/statsmodels/tools/_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
  import pandas.util.testing as tm
In [2]:
import plotly.graph_objs as go
In [3]:
import os
fpath = os.path.join(os.path.dirname(pyActigraphy.__file__),'tests/data/')
In [4]:
raw = pyActigraphy.io.read_raw_rpx(fpath+'test_sample.csv', period='6 days', language='FR')

Epoch-by-epoch rest/activity scoring algorithms

All based on the same idea;

  • use a rolling window over the data and convolute them with a “gaussian”-like kernel.

  • if the resulting data is above a predefined threshold, classify as activity. Rest otherwise…

My opinion regarding these algorithms

The weights of these “gaussian”-like kernels as well as the applied thresholds are very much device and population dependant. Using these algorithms would require to adapt these parameters to the population under study… Possible but seldom done.

In [5]:
layout = go.Layout(title="Rest/Activity detection",xaxis=dict(title="Date time"), yaxis=dict(title="Counts/period"), showlegend=False)

Cole-Kripke

Disclaimer: in the original paper, they define “rescoring” rules… Not implemented in the package.

In [6]:
help(raw.CK)
Help on method CK in module pyActigraphy.sleep.scoring_base:

CK(scale=1e-05, window=array([ 400,  600,  300,  400, 1400,  500,  350,    0,    0], dtype=int32), threshold=1.0) method of pyActigraphy.io.rpx.rpx.RawRPX instance
    Cole&Kripke algorithm for sleep-wake identification.

    Algorithm for automatic sleep scoring based on wrist activity,
    developped by Cole, Kripke et al [1]_.


    Parameters
    ----------
    scale: float, optional
        Scale parameter P.
        Default is 0.00001.
    window: np.array
        Array of weighting factors, :math:`W_{i}`.
        Default is [400, 600, 300, 400, 1400, 500, 350, 0, 0]
    threshold: float, optional
        Threshold value for scoring sleep/wake.
        Default is 1.0.

    Returns
    -------
    ck: pandas.core.Series
        Time series containing the `D` scores (0: sleep, 1: wake) for each
        epoch.

    Notes
    -----

    The output variable D of the CK algorithm is defined in [1]_ as:

    .. math::

        D = P*(
            [W_{-4},\dots,W_{0},\dots,W_{+2}]
            \cdot
            [A_{-4},\dots,A_{0},\dots,A_{+2}])

    with:

    * D < 1 == sleep, D >= 1 == wake;
    * P, scale factor;
    * :math:`W_{0},W_{-1},W_{+1},\dots`, weighting factors for the present
      minute, the previous minute, the following minute, etc.;
    * :math:`A_{0},A_{-1},A_{+1},\dots`, activity scores for the present
      minute, the previous minute, the following minute, etc.


    References
    ----------

    .. [1] Cole, R. J., Kripke, D. F., Gruen, W., Mullaney, D. J.,
           & Gillin, J. C. (1992). Automatic Sleep/Wake Identification
           From Wrist Activity. Sleep, 15(5), 461–469.
           http://doi.org/10.1093/sleep/15.5.461

    Examples
    --------

In [7]:
CK = raw.CK()
In [8]:
layout.update(yaxis2=dict(title='Classification',overlaying='y',side='right'), showlegend=True);
In [9]:
go.Figure(data=[
    go.Scatter(x=raw.data.index.astype(str),y=raw.data, name='Data'),
    go.Scatter(x=CK.index.astype(str),y=CK, yaxis='y2', name='CK')
], layout=layout)

Sadeh’s and Scripps’ algorithms

In [10]:
sadeh = raw.Sadeh()
scripps = raw.Scripps()
In [11]:
go.Figure(data=[
    go.Scatter(x=raw.data.index.astype(str),y=raw.data, name='Data'),
    go.Scatter(x=sadeh.index.astype(str),y=sadeh, yaxis='y2', name='Sadeh'),
    go.Scatter(x=scripps.index.astype(str),y=scripps, yaxis='y2', name='Scripps')
], layout=layout)

Oakley’s algorithm

This is the sleep/wake scoring algorithm used by the Actiware software from Respironics. The various activity thresholds have thus been tailored for Actiwatch devices. Be careful when applying this algorithm on data acquired with other devices.

In [12]:
help(raw.Oakley)
Help on method Oakley in module pyActigraphy.sleep.scoring_base:

Oakley(threshold=40) method of pyActigraphy.io.rpx.rpx.RawRPX instance
    Oakley's algorithm for sleep/wake scoring.

    Algorithm for automatic sleep/wake scoring based on wrist activity,
    developed by Oakley [1]_.


    Parameters
    ----------
    threshold: float or str, optional
        Threshold value for scoring sleep/wake. Can be set to "automatic"
        (cf. Notes).
        Default is 40.

    Returns
    -------

    oakley: pandas.core.Series
        Time series containing scores (1: sleep, 0: wake) for each
        epoch.


    Notes
    -----

    The output variable O of Oakley's algorithm is defined as:

    .. math::

        O = (
            [W_{-2},W_{-1},W_{0},W_{+1},W_{+2}]
            \cdot
            [A_{-2},A_{-1},A_{0},A_{+1},A_{+2}])

    with:

    * O <= threshold == sleep, 0 > threshold == wake;
    * :math:`W_{0},W_{-1},W_{+1},\dots`, weighting factors for the present
      epoch, the previous epoch, the following epoch, etc.;
    * :math:`A_{0},A_{-1},A_{+1},\dots`, activity scores for the present
      epoch, the previous epoch, the following epoch, etc.


    The current implementation of this algorithm follows the description
    provided in the instruction manual of the Actiwatch Communication and
    Sleep Analysis Software (Respironics, Inc.) [2]_ :

    * 15-second sampling frequency:
        .. math::

            W_{-8} &= \ldots = W_{-5} = W_{+5} = \ldots = W_{+8} = 1/25 \\
            W_{-4} &= \ldots = W_{-1} = W_{+1} = \ldots = W_{+4} = 1/5 \\
            W_{0} &= 4

    * 30-second sampling frequency:
        .. math::

            W_{-4} &= W_{-3} = W_{+3} = W_{+4} = 1/25 \\
            W_{-2} &= W_{-1} = W_{+1} = W_{+2} = 1/5 \\
            W_{0} &= 2

    * 60-second sampling frequency:
        .. math::

            W_{-2} &= W_{+2} = 1/25 \\
            W_{-1} &= W_{+1} = 1/5 \\
            W_{0} &= 1

    * 120-second sampling frequency:
        .. math::

            W_{-1} &= W_{+1} = 1/8 \\
            W_{0} &= 1/2


    The *Automatic Wake Threshold Value* calculation is this [2]_:

    1. Sum the activity counts for all epochs of the data set.
    2. Count the number of epochs scored as MOBILE for the data set
       (the definition of MOBILE follows).
    3. Compute the MOBILE TIME (number of epochs scored as MOBILE from
       step 2 multiplied by the Epoch Length) in minutes.
    4. Compute the Auto Wake Threshold = ((sum of activity counts from
       step 1) divided by (MOBILE TIME from step 3)) multiplied by 0.88888.


    *Definition of Mobile* [2]_:

    An epoch is scored as MOBILE if the number of activity counts recorded
    in that epoch is greater than or equal to the epoch length in 15-second
    intervals. For example,there are four 15-second intervals for a
    1-minute epoch length; hence, the activity value in an epoch must be
    greater than, or equal to, four, to be scored as MOBILE.


    References
    ----------

    .. [1] Oakley, N.R. Validation with Polysomnography of the Sleepwatch
           Sleep/Wake Scoring Algorithm Used by the Actiwatch Activity
           Monitoring System; Technical Report; Mini-Mitter: Bend, OR, USA,
           1997
    .. [2] Instruction manual, Actiwatch Communication and Sleep Analysis
           Software
           (https://fccid.io/JIAAWR1/Users-Manual/USERS-MANUAL-1-920937)

Medium threshold

A threshold of 40 corresponds to a medium wake threshold.

In [13]:
oakley = raw.Oakley(threshold=40)

Automatic threshold

It is also possible to compute a threshold value automatically, based on activity data.

In [14]:
oakley_auto = raw.Oakley(threshold='automatic')
In [15]:
go.Figure(data=[
    go.Scatter(x=raw.data.index.astype(str),y=raw.data, name='Data'),
    go.Scatter(x=oakley.index.astype(str),y=sadeh, yaxis='y2', name='Oakley (thr: medium)'),
    go.Scatter(x=oakley_auto.index.astype(str),y=scripps, yaxis='y2', name='Oakley (thr: automatic)')
], layout=layout)

Consolidated activity/rest period detection

Crespo’s algorithm

This is a threshold-based algorithm that used morphological filters to “clean” short periods of activity (rest) surrounded by periods of rest (acitivity).

In [16]:
help(raw.Crespo)
Help on method Crespo in module pyActigraphy.sleep.scoring_base:

Crespo(zeta=15, zeta_r=30, zeta_a=2, t=0.33, alpha='8h', beta='1h', estimate_zeta=False, seq_length_max=100, verbose=False) method of pyActigraphy.io.rpx.rpx.RawRPX instance
    Crespo algorithm for activity/rest identification

    Algorithm for automatic identification of activity-rest periods based
    on actigraphy, developped by Crespo et al. [1]_.

    Parameters
    ----------
    zeta: int, optional
        Maximum number of consecutive zeroes considered valid.
        Default is 15.
    zeta_r: int, optional
        Maximum number of consecutive zeroes considered valid (rest).
        Default is 30.
    zeta_a: int, optional
        Maximum number of consecutive zeroes considered valid (active).
        Default is 2.
    t: float, optional
        Percentile for invalid zeroes.
        Default is 0.33.
    alpha: str, optional
        Average hours of sleep per night.
        Default is '8h'.
    beta: str, optional
        Length of the padding sequence used during the processing.
        Default is '1h'.
    estimate_zeta: bool, optional
        If set to True, zeta values are estimated from the distribution of
        ratios of the number of series of consecutive zeroes to
        the number of series randomly chosen from the actigraphy data.
        Default is False.
    seq_length_max: int, optional
        Maximal length of the aforementioned random series.
        Default is 100.
    verbose: bool, optional
        If set to True, print the estimated values of zeta.
        Default is False.

    Returns
    -------
    crespo : pandas.core.Series
        Time series containing the estimated periods of rest (0) and
        activity (1).

    References
    ----------

    .. [1] Crespo, C., Aboy, M., Fernández, J. R., & Mojón, A. (2012).
           Automatic identification of activity–rest periods based on
           actigraphy. Medical & Biological Engineering & Computing, 50(4),
           329–340. http://doi.org/10.1007/s11517-012-0875-y

    Examples
    --------

        >>> import pyActigraphy
        >>> rawAWD = pyActigraphy.io.read_raw_awd(fpath + 'SUBJECT_01.AWD')
        >>> crespo = rawAWD.Crespo()
        >>> crespo
        2018-03-26 14:16:00    1
        2018-03-26 14:17:00    0
        2018-03-26 14:18:00    0
        (...)
        2018-04-06 08:22:00    0
        2018-04-06 08:23:00    0
        2018-04-06 08:24:00    1
        Length: 15489, dtype: int64

In [17]:
crespo = raw.Crespo()
In [18]:
crespo_6h = raw.Crespo(alpha='6h')
In [19]:
crespo_zeta = raw.Crespo(estimate_zeta=True)
In [20]:
go.Figure(data=[
    go.Scatter(x=raw.data.index.astype(str),y=raw.data, name='Data'),
    go.Scatter(x=crespo.index.astype(str),y=crespo, yaxis='y2', name='Crespo'),
    go.Scatter(x=crespo_6h.index.astype(str),y=crespo_6h, yaxis='y2', name='Crespo (6h)'),
    go.Scatter(x=crespo_zeta.index.astype(str),y=crespo_zeta, yaxis='y2', name='Crespo (Automatic)')
], layout=layout)

Like for the other algorithms, the default parameters correspond to those described in the original papers, which, most likely, have been tuned to the population used to validate the algorithm.

Since the output is made of consolidated periods, it is possible to return the offset and onset times of each period:

In [21]:
aot = raw.Crespo_AoT()
In [22]:
aot
Out[22]:
(DatetimeIndex(['2015-02-04 12:45:00', '2015-02-05 06:13:00',
                '2015-02-06 07:44:00', '2015-02-07 05:45:00',
                '2015-02-08 06:22:00', '2015-02-09 07:09:00',
                '2015-02-10 05:56:00', '2015-02-10 11:45:00'],
               dtype='datetime64[ns]', freq=None),
 DatetimeIndex(['2015-02-04 11:46:00', '2015-02-04 22:07:00',
                '2015-02-05 19:49:00', '2015-02-06 19:43:00',
                '2015-02-07 21:28:00', '2015-02-08 19:27:00',
                '2015-02-09 19:42:00', '2015-02-10 08:44:00'],
               dtype='datetime64[ns]', freq=None))
In [23]:
aot[0]-aot[1]
Out[23]:
TimedeltaIndex(['00:59:00', '08:06:00', '11:55:00', '10:02:00', '08:54:00',
                '11:42:00', '10:14:00', '03:01:00'],
               dtype='timedelta64[ns]', freq=None)

Roenneberg’s algorithm

Also threshold-based but uses correlations with test series of various lengths to find the consolidated period that best matches best the data.

In [24]:
help(raw.Roenneberg)
Help on method Roenneberg in module pyActigraphy.sleep.scoring_base:

Roenneberg(trend_period='24h', min_trend_period='12h', threshold=0.15, min_seed_period='30Min', max_test_period='12h', n_succ=3) method of pyActigraphy.io.rpx.rpx.RawRPX instance
    Automatic sleep detection.

    Identification of consolidated sleep episodes using the
    algorithm developped by Roenneberg et al. [1]_.

    Parameters
    ----------
    trend_period: str, optional
        Time period of the rolling window used to extract the data trend.
        Default is '24h'.
    min_trend_period: str, optional
        Minimum time period required for the rolling window to produce a
        value. Values default to NaN otherwise.
        Default is '12h'.
    threshold: float, optional
        Fraction of the trend to use as a threshold for sleep/wake
        categorization.
        Default is '0.15'
    min_seed_period: str, optional
        Minimum time period required to identify a potential sleep onset.
        Default is '30Min'.
    max_test_period : str, optional
        Maximal period of the test series.
        Default is '12h'
    n_succ : int, optional
        Number of successive elements to consider when searching for the
        maximum correlation peak.

    Returns
    -------
    rbg : pandas.core.Series
        Time series containing the estimated periods of rest (1) and
        activity (0).

    References
    ----------

    .. [1] Roenneberg, T., Keller, L. K., Fischer, D., Matera, J. L.,
           Vetter, C., & Winnebeck, E. C. (2015). Human Activity and Rest
           In Situ. In Methods in Enzymology (Vol. 552, pp. 257-283).
           http://doi.org/10.1016/bs.mie.2014.11.028

    Examples
    --------

In [25]:
roenneberg = raw.Roenneberg()
roenneberg_thr = raw.Roenneberg(threshold=0.25, min_seed_period='15min')
/Users/ghammad/Work/CRC/pyActigraphy/pyActigraphy/sleep/scoring/utils.py:99: RuntimeWarning:

invalid value encountered in greater

In [26]:
go.Figure(data=[
    go.Scatter(x=raw.data.index.astype(str),y=raw.data, name='Data'),
    go.Scatter(x=roenneberg.index.astype(str),y=roenneberg, yaxis='y2', name='Roenneberg'),
    go.Scatter(x=roenneberg_thr.index.astype(str),y=roenneberg_thr, yaxis='y2', name='Roenneberg (Thr:0.25)')
], layout=layout)
In [27]:
aot = raw.Roenneberg_AoT(threshold=0.25, min_seed_period='15min')
/Users/ghammad/Work/CRC/pyActigraphy/pyActigraphy/sleep/scoring/utils.py:99: RuntimeWarning:

invalid value encountered in greater

In [28]:
aot
Out[28]:
(DatetimeIndex(['2015-02-04 13:03:00', '2015-02-05 08:16:00',
                '2015-02-05 10:56:00', '2015-02-05 15:08:00',
                '2015-02-05 15:43:00', '2015-02-06 03:31:00',
                '2015-02-06 07:36:00', '2015-02-06 11:15:00',
                '2015-02-06 19:32:00', '2015-02-06 20:06:00',
                '2015-02-07 07:47:00', '2015-02-07 11:01:00',
                '2015-02-07 13:37:00', '2015-02-07 18:58:00',
                '2015-02-08 07:35:00', '2015-02-08 14:31:00',
                '2015-02-09 07:20:00', '2015-02-09 11:06:00',
                '2015-02-10 06:39:00', '2015-02-10 07:32:00'],
               dtype='datetime64[ns]', name='Date_Time', freq=None),
 DatetimeIndex(['2015-02-04 12:45:00', '2015-02-04 21:07:00',
                '2015-02-05 10:32:00', '2015-02-05 14:41:00',
                '2015-02-05 15:24:00', '2015-02-05 19:25:00',
                '2015-02-06 03:50:00', '2015-02-06 10:56:00',
                '2015-02-06 18:58:00', '2015-02-06 19:44:00',
                '2015-02-06 20:37:00', '2015-02-07 10:42:00',
                '2015-02-07 13:19:00', '2015-02-07 18:36:00',
                '2015-02-07 21:00:00', '2015-02-08 14:01:00',
                '2015-02-08 20:10:00', '2015-02-09 10:42:00',
                '2015-02-09 18:55:00', '2015-02-10 06:48:00'],
               dtype='datetime64[ns]', name='Date_Time', freq=None))
In [29]:
aot[0]-aot[1]
Out[29]:
TimedeltaIndex(['00:18:00', '11:09:00', '00:24:00', '00:27:00', '00:19:00',
                '08:06:00', '03:46:00', '00:19:00', '00:34:00', '00:22:00',
                '11:10:00', '00:19:00', '00:18:00', '00:22:00', '10:35:00',
                '00:30:00', '11:10:00', '00:24:00', '11:44:00', '00:44:00'],
               dtype='timedelta64[ns]', name='Date_Time', freq=None)

Et voilà! Easy, isn’t it?